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The structural, dynamical, and thermodynamical properties of diamond, graphite and layered 
derivatives (graphene, rhombohedral graphite) are computed using a combination of density- 
functional theory (DFT) total-energy calculations and density-functional perturbation theory 
(DFPT) lattice dynamics at the GGA-PBE level. Overall, very good agreement is found for the 
structural properties and phonon dispersions, with the exception of the c/a ratio in graphite and 
the associated elastic constants and phonon dispersions. Both the C33 elastic constant and the F to 
A phonon dispersions are brought to close agreement with available data once the experimental c/a 
is chosen for the calculations. The thermal expansion, the temperature dependence of the elastic 
moduli and the specific heat have been calculated via the quasi-harmonic approximation. Graphite 
shows a distinctive in-plane negative thermal-expansion coefficient that reaches the minimum around 
room temperature, in very good agreement with experiments. Thermal contraction in graphene is 
found to be three times as large; in both cases, ZA acoustic modes are shown to be responsible for 
the contraction, in a direct manifestation of the membrane effect predicted by Lifshitz over fifty 
years ago. 
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I. INTRODUCTION 

The extraordinary variety of carbon allotropes, as well 
as their present and potential applications in such di- 
verse fields as nanoelectronicai or bioengineeringl, gives 
them a special place among all the elements. Even ex- 
cluding fuUerenes, nanotubes, and their derivatives, sin- 
gle crystalline diamond, graphite and graphene (i.e. a 
single graphite layer) still lack a complete characteri- 
zation of their thermodynamic sta bilit y under a broad 
range of conditions (see e.g. Refs. |3l33al3 and citations 
therein). In this respect, vibrational properties play a 
crucial role in determining the thermodynamic proper- 
ties of the bulk. Indeed, diamond being a wide band 
gap material (Eg= 5.5 eV), electronic excitations do not 
account for thermal properties up to high temperatures. 
Graphite and graphene are semi-metals, but the gap van- 
ishes only at the K point where the two massless bands 
cross (see e.g. Ref. @); thus, electronic excitations can 
also be neglected in these materials, and the phonon dis- 
persions provide all the information that is needed to cal- 
culate thermodynamical quantities such as the thermal 
expansion or specific heat. 

The aim of this paper is to provide a con- 
verged, accurate determination of the structural, dy- 
namical, and thermodynamical properties of diamond, 
graphite, graphene and rhombohedral graphite from first- 
principles calculations. Although the phonon spectrum 
of diamond and its thermal properties have been stud- 
ied extensively with experimentsSiifl and calculationsii, 
the phonon spectrum of graphite is still under ac- 
tive investigation^^'"'^'^, as well as its thermal proper- 
ties. Graphite in-plane thermal expansion has long been 
recognized to be negativ o^^i^^ , and it has even been 



suggestedSii^ that this may be due to the internal stresses 
related to the large expansion in the c direction (Poisson 
effect). To resolve some of the open questions, and to 
provide a coherent picture for these materials, we used 
extensive ab-initio density-functional theory (DFT) and 
density-functional perturbation theory (DFPT)iS4ii cal- 
culations. DFT is a very efficient and accurate tool to 
obtain ground-state and linear-response properties, es- 
pecially when paired with plane-wave basis sets, which 
easily allow to reach full convergence with respect to 
basis size, and ultrasoft pseudo-potentials^ for optimal 
performance and transferability. We adopted the PBE- 
GGA^^ exchange-correlation functional, at variance with 
most of the early studies on diamondiiiSS*2i and espe- 
cially graphitei2*2242iSii2^i2&, which have been performed 
using the local density approximation (LDA). GGA cal- 
culations have appeared_ mostly for the cases of dia- 
mond (GGA-PBE, Ref. 21) and graphene (GGA-PBE, 
Refs. tl^ilSi ), wit h some data for graphite appearing in 



Refs.liiailili (GGA-PBE). DEPTiSJi is then used to 
compute the phonon frequencies at any arbitrary wave- 
vector, without having to resort to the use of supercells. 
The vibrational free energy is calculated in the quasi- 
harmonic approximation (QIIA)ii^, to predict finite- 
temperature lattice properties such as thermal expansion 
and specific heat. 

To the best of our knowledge, this is the first study on 
the thermal properties of graphite or graphene from first- 
principles. For the case of diamond and graphene, calcu- 
lations are fully ab-initio and do not require any experi- 
mental input. For the case of graphite and rhombohedral 
graphite we argue that the use of the experimental c/a 
greatly improves the agreement with experimental data. 
This experimental input is required since DFT, in its 



2 



current state of development, yields poor predictions for 
the interlayer interactions, dominated by Van Der Waals 
dispersion forces not well described by local or semi-local 
exchange correlation functionals (see Refs. "sT and's? for 
details; the agreement between LDA predictions and ex- 
perimental results for the cja ratio is fortuitous). It is 
found that the weak interlayer bonding has a small influ- 
ence on most of the properties studied and that forcing 
the experimental cja corrects almost all the remaining 
ones. This allowed us to obtain results for all the materi- 
als considered that are in very good agreement with the 
available experimental data. 

The article is structured as follows. We give a brief 
summary of our approach and definitions and introduce 
DFPT and the QHA in Section HTJ Our ground-state, 
zero-temperature results for diamond, graphite, graphene 
and rhombohedral graphite are presented in Section IlIII 
Lattice parameters and elastic constants from the equa- 
tions of state in subsection llll Al phonon frequencies and 
vibrational density of states in siibsection llll Bl and first- 
principles, linear-response interatomic force constants in 
subsection IIII CI The lattice thermal properties, such 
as thermal expansion, mode Griineisen parameters, and 
specific heat as obtained from the vibrational free energy 
are presented in section HVI Section M contains our final 
remarks. 
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Here R (R') is a Bravais lattice vector, i (j) indicates 
the i*'' (j*'') atom of the unit cell, and a{/3) repre- 
sents the cartesian components. C^°"^j are the second 
derivatives^^ of Ewald sums corresponding to the ion- 
ion repulsion potential, while the electronic contributions 
^ai'pj second derivatives of the electron-electron 

and electron-ion terms in the ground state energy. From 
the Hellmann-FeynmaniZ. theorem one obtains: 
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II. THEORETICAL FRAMEWORK 

A. Basics of Density-Functional Perturbation 
Theory 

In density-functional thcoryS*^ the ground state elec- 
tronic density and wavefunctions of a crystal are found by 
solving self-consistently a set of one-electron equations. 
In atomic units (used throughout the article), these are 

[-UJ^ + VsCF{vm^)=e^\i^^). (la) 



^-.(r)^/^^V + ^ + y.„(r), (lb) 
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where /(sp—Si) is the occupation function, ep the Fermi 
energy, E^c the exchange-correlation functional (approx- 
imated by GGA-PBE in our case), n{r) the electronic- 
density, and Vio„(r) the ionic core potential (actually a 
sum over an array of pseudo-potentials). 

Once the unperturbed ground state is determined, 
phonon frequencies can be obtained from the interatomic 
force constants, i.e. the second derivatives at equilibrium 



(where the dependence of both n(r) and Vion(r) on the 
displacements has been omitted for clarity, and T^on(r) 
is considered local). 

It is seen that the electronic contribution can be ob- 
tained from the knowledge of the linear response of the 
system to a displacement. The key assumption is then 
the Born-Oppenheimer approximation which views a lat- 
tice vibration as a static perturbation on the electrons. 
This is equivalent to say that the response time of the 
electrons is much shorter than that of ions, that is, each 
time ions are slightly displaced by a phonon, electrons 
instantaneously rearrange themselves in the state of min- 
imum energy of the new ionic configuration. Therefore, 
static linear response theory can be applied to describe 
the behavior of electrons upon a vibrational excitation. 

For phonon calculations, we consider a periodic pertur- 
bation AVion of wave-vector q, which modifies the self- 
consistent potential Vscf by an amount AVscf- The 
linear response in the charge density An(r) can be found 
using first-order perturbation theory. If we consider its 
Fourier transform An(q + G), and calling V'o.k the one- 
particle wavefunction of an electron in the occupied band 
"o" at the point k of the Brillouin zone (and £o,k the cor- 
responding eigenvalue), one can get a self-consistent set 
of linear equations similar to Eqs. Hlallbllc^ .2^: 



(£o,k + 2 V - ■^5CF(r))AVo,k+q = P^+'^AV^cF^oM 
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An(q + G) = - ^(V^e,k|e-*('>+°)--A''+'^|A^„.k+q) 
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of independent harmonic oscillators. In a straightforward 
manner, it can be shown'^'' that: 
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pk+q j-gfgj-g the projector on the empty-state manifold 
at k + q, V to the total crystal volume, and G to any 
reciprocal lattice vector. Note that the linear response 
contains only Fourier components of wave vector q + G, 
so we added a superscript q to AVg^jp. We have implic- 
itly assumed for simplicity that the crystal has a band 
gap and that pseudo-potentials are local, but the general- 
ization to metals?^ and to non-local pseudo-potentials^ - 
are all well established (see Ref. for a detailed and 
complete review of DFPT). 

Linear-response theory allows us to calculate the re- 
sponse to any periodic perturbation; i.e. it allows direct 
access to the dynamical matrix related to the interatomic 
force constants via a Fourier transform: 



« (q) = -j^== E « (R) e~^^ ^ (5) 

(where Mi is the mass of the i*'' atom). 

Phonon frequencies at any q are the solutions of the 
eigenvalue problem: 
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In practice, one calculates the dynamical matrix on a rel- 
atively coarse grid in the Brillouin zone (say, a 8 x 8 x 8 
grid for diamond), and obtains the corresponding inter- 
atomic force constants by inverse Fourier transform (in 
this example it would correspond to a 8 x 8 x 8 super- 
cell in real space). Finally, the dynamical matrix (and 
phonon frequencies) at any q point can be obtained by 
Fourier interpolation of the real-space interatomic force 
constants. 



B. Thermodynamical properties 

When no external pressure is applied to a crystal, the 
equilibrium structure at any temperature T can be found 
by minimizing the Helmholtz free energy F{{ai},T) = 
U — TS with respect to all its geometrical degrees of free- 
dom {ui}. If now the crystal is supposed to be perfectly 
harmonic, F is the sum of the ground state total energy 
and the vibrational free energy coming from the parti- 
tion function (in the canonical ensemble) of a collection 
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where E{{ai}) is the ground state energy and the sums 
run over all the Brillouin zone wave- vectors and the band 
index j of the phonon dispersion. The second term in the 
right hand side of Eq. O is the zero-point motion. 

If anharmonic effects are neglected, the phonon fre- 
quencies do not depend on lattice parameters, therefore 
the free energy dependence on structure is entirely con- 
tained in the ground state equation of state E{{ai\). 
Consequently the structure does not depend on temper- 
ature in a harmonic crystal. 

Thermal expansion is recovered by introducing in 
Eq. Q the dependence of the phonon frequencies on the 
structural parameters {a^}; direct minimization of the 
free energy 
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provides the equilibrium structure at any temperature T. 
This approach goes under the name quasi-harmonic ap- 
proximation (QHA) and has been applied successfully to 
many bulk systemsii*^2i^. The linear thermal expansion 
coefficients of the cell dimensions of a lattice arc then 
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The Griineisen formalism^^" assumes a linear dependence 
of the phonon frequencies on the three orthogonal cell 
dimensions {a^}; developing the ground state energy up 
to second order, (thanks to the equation of state at T = 

Oif), one can get from the condition = the 

alternative expression 
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We follow here the formalism of Ref. [IJ c„(q, j) is the 
contribution to the specific heat from the mode (q, j). 
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Sik is the elastic compliance matrix, and the subscript 
"0" indicates a quantity taken at the ground state lattice 
parameter. The Griineisen parameter of the mode (q, j) 
is by definition 



7fc(q,i) = 



^o,q,i dak 
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For a structure which depends only on one lattice param- 
eter a (e.g. diamond or graphene) one then gets for the 
linear thermal expansion coefficient 
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where Bq is defined by Bq = Vb§^ {V represents the 
volume of a three-dimensional crystal such as diamond 
or the surface of a two-dimensional one like graphene) , d 
is the number of dimensions (d = 3 for diamond, d — 2 
for graphene), and Vq is the volume (or the surface) at 
equilibrium. 

In the case of graphite there are two lattice parameters: 
a in the basal plane and c perpendicular to the basal 
plane, so that one gets 



2w, 



da 



513- 



dc 



(13a) 



wo,qj da 



-Co duq^. 



dc 



(13b) 



The mode Griineisen parameters provide useful insight 
to the thermal expansion mechanisms. They are usu- 
ally positive, since phonon frequencies decrease when the 
solid expands, although some negative mode Griineisen 
parameters for low-frequency acoustic modes can arise 
and sometimes compete with the positive ones, giving a 
negative thermal expansion at low temperatures, when 
only the lowest acoustic modes can be excited. 

Finally, the heat capacity of the unit cell at constant 
volume can be obtained from C„ 
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C. Computational details 

All the calculations that follow were performed us- 
ing the ESPRESSO^^ package, which is a fuU ab- 
initio DFT and DEFT code available under the 
GNU Fublic License^. We used a plane-wave ba- 
sis set, ultrasoft pseudo-potentials''^ from the standard 
distribution*'* (generated using a modified RRKJ^S. ap- 
proach), and the generalized gradient approximation 
(GGA) for the exchange-correlation functional in its PBE 
parameterization*^. We also used the local density ap- 
proximation (LDA) in order to compare some results be- 
tween the two functionals. In this case the parameteriza- 
tion used was the one proposed by Ferdew and Zunger— . 

For the semi-metallic graphite and graphene cases, we 
used 0.03 Ryd of cold smearing^. We carefully and ex- 
tensively checked the convergence in the energy differ- 
ences between different configurations and the phonon 
frequencies with respect to the wavefunction cutoff, 
the dual (i.e. the ratio between charge density cut- 
off and wavefunction cutoff), the k-point sampling of 
the Brillouin zone, and the interlayer vacuum spacing 
for graphene. Energy differences were converged within 
5 meV/atom or better, and phonon frequencies within 
1 — 2 cm~*. In the case of graphite and graphene phonon 
frequencies were converged with respect to the k-point 
sampling after having set the smearing parameter at 0.03 
Ryd. Besides, values of the smearing between 0.02 Ryd 
and 0.04 Ryd did not change the frequencies by more 
than 1 — 2 cm~*. 

In a solid, translational invariance guaranties that 
three phonon frequencies at F will go to zero. In our 
GGA-FBE DEFT formalism this condition is exactly sat- 
isfied only in the limit of infinite k-point sampling and 
full convergence with the plane-wave cutoff. For the case 
of graphene and graphite we found in particular that an 
exceedingly large cutoff (100 Ryd) and dual (28) would be 
needed to recover phonon dispersions (especially around 
r and the F — ^ branch) with the tolerances mentioned; 
on the other hand, application of the acoustic sum rule 
(i.e. forcing the translational symmetry on the inter- 
atomic force constants) allows us to recover these highly 
converged calculations above with a more reasonable cut- 
off and dual. 

Finally, the cutoffs we used were 40 Ryd for the wave- 
functions in all the carbon materials presented, with du- 
als of 8 for diamond and 12 for graphite and graphene. 
We used a 8 x 8 x 8 Monkhorst-Pack k-mesh for diamond, 
16 X 16 X 8 for graphite, 16 x 16 x 4 for rhombohedral 
graphite and 16 x 16 x 1 for graphene. All these meshes 
were unshifted (i.e. they do include F). Dynamical ma- 
trices were initially calculated on a 8 x 8 x 8 q-points 
mesh for diamond, 8x8x4 for graphite, 8 x 8 x 2 for 
rhombohedral graphite and 16 x 16 x 1 for graphene. 

Finally, integrations over the Brillouin zone for the vi- 
brational free energy or the heat capacity were done us- 
ing phonon frequencies that were Fourier interpolated on 
much finer meshes. The phonon frequencies were usually 
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TABLE I: Equilibrium lattice parameter ao and bulk mod- 
ulus Bo of diamond at the ground state (GS) and at 300 K 
(see Section HVt . compared to experimental values. 



Present calculation Experiment (300 K) 


Lattice constant ao 


6.743 (GS) 6.740 " 


(a.u.) 


6.769 (300 K) 


Bulk modulus Bo 


432 (GS) 442 ± 2 " 


(GPa) 


422 (300 K) 


"Ref.ra 




''Ref.lig 





computed at several lattice parameters and the results 
interpolated to get their dependence on lattice constants. 

A final remark is that we were careful to use the same 
parameters (cutoffs, k-points sampling, smearing, etc.) 
in the determination of the ground state equation of state 
and that of the phonon frequencies, since these two terms 
need to be added in the free energy expression. 



III. ZERO-TEMPERATURE RESULTS 
A. Structural and elastic properties 

We performed ground state total-energy calculations 
on diamond, graphite, and graphene over a broad range 
of lattice parameters. The potential energy surface can 
then be fitted by an appropriate equation of state. The 
minimum gives the ground state equihbrium lattice pa- 
rameter(s). The second derivatives at that minimum are 
related to the bulk modulus or elastic constants. 

For the case of diamond we chose the Birch equation 
of stated (up to the fourth order) to fit the total energy 
vs. the lattice constant a: 
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where Bq is the bulk modulus, Vq the primitive cell vol- 
ume (Vb = ^ here) and A and B are fit parameters. 
The Murnaghan equation of state or even a polynomial 
would fit equally well the calculations around the min- 
imum of the curve. A best fit of this equation on our 
data gives us both the equilibrium lattice parameter and 
the bulk modulus; our results are summarized in Tabled 
The agreement with the experimental values is very good, 
even after the zero-point motion and thermal expansion 
are added to our theoretical predictions (see Section llV)l . 

The ground state equation of state of graphene was fit- 
ted by a 4'^ order polynomial, and the minimum found 
for a = 4.654 a.u., which is very close to the experimen- 
tal in-plane lattice parameter of graphite. The graphite 
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FIG. 1: Contour plot of the ground state energy of graphite 
as a function of a and c/a (isoenergy contours are not equidis- 
tant). 
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FIG. 2: Ground state energy of graphite as a function of 
c/a at fixed a — 4.65 a.u.. The theoretical (PBE) and the 
experimental c/a axe shown. The zero of energy has been set 
to the PBE minimum. 



equation of state was fitted by a two-dimensional 4*^ or- 
der polynomial of variables a and c. To illustrate the 
very small dependence of the ground state energy with 
the c/a ratio, we have plotted the results of our calcu- 
lations over a broad range of lattice constants in Figs. 
and 121 

A few elastic constants can be obtained from the sec- 
ond derivatives of this energy^: 
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Tetragonal shear modulus C' = - [(Cn + C12 



Bulk modulus B, 







+ 2C33 - 4Ci3] (16b) 
C33(C'ii + C12) — 26*^3 



6C* 



(16c) 



We summarize all our LDA and GGA results in Ta- 
ble For LDA, both the lattice parameter oq and the 
Co/ao ratio are very close to experimental data. Elastic 
constants were calculated fully from first-principles, in 
the sense that the second derivatives of the energy were 
taken at the theoretical LDA ao and co, and that only 
these theoretical values were used in Eqs. (|16a|) . Elastic 
constants are found in good agreement with experiments, 
except for the case of C13 which comes out as negative 
(meaning that the Poisson's coefhcicnt would be nega- 
tive). 

Fully theoretical GGA results (second column of Ta- 
ble ^ compare poorly to experimental data except for 
the ao lattice constant, in very good agreement with 
experiments. Using the experimental value for cq in 
Eqs. ((TB|) improves only the value of Cn -I- C12 (third 
column of Table ^ . Most of the remaining disagree- 
ment is related to the poor value obtained for c/a; if the 
second derivatives in Eqs. H16a|l are taken at the exper- 
imental value for c/a all elastic constants are accurately 
recovered except for C13 (fourth column of Table ITU) . 

In both LDA and GGA, errors arise from the fact that 
Van Der Waals interactions between graphitic layers are 
poorly described. These issues can still be addressed 
within the framework of DFT (as shown by Langreth and 
collaborators, Ref. Isil) at the cost of having a non-local 
exchange-correlation potential. 

Zero-point motion and finite-temperature effects will 
be discussed in detail in Section Hvl 



B. Phonon dispersion curves 

We have calculated the phonon dispersion relations 
for diamond, graphite, rhombohedral graphite and 
graphene. For diamond and graphene, we used the the- 
oretical lattice parameter. For graphite, we either used 
the theoretical c/a or the experimental one (c/a = 2.725). 
We will comment extensively in the following on the role 
of c/a on our calculated properties. 

Finally we also calculated the phonon dispersions for 
rhombohedral graphite, which differs from graphite only 
in the stacking of the parallel layers: in graphite the 
stacking is ABABAB while it is ABCABC in rhombohe- 
dral graphite, and the latter unit cell contains six atoms 
instead of four. We therefore used the same in-plane lat- 
tice parameter and same interlayer distance as in graphite 
(that is, a ^ ratio multiplied by | ). Results are presented 
in Figs.|3121[Sliniand[7| together with the experimental 



data. 

In Table UTTl and Hvl we summarize our results at high- 
symmetry points and compare them with experimental 
data. In diamond, GGA produces softer modes than 
LDAii on the whole (as expected) , particularly at F (op- 
tical mode) and in the optical F-X branches. For these, 
the agreement is somehow better in LDA; on the other 
hand the whole F-L dispersion is overestimated by LDA. 

The results on graphite require some comments. In 
Table IIVI and Figs. 01 El El and modes are classified 
as follow: L stands for longitudinal polarization, T for 
in-plane transversal polarization and Z for out-of-plane 
transversal polarization. For graphite, a prime (as in 
LO') indicates an optical mode where the two atoms in 
each layer of the unit cell oscillate together and in phase 
opposition to the two atoms of the other layer. A non- 
primed optical mode is instead a mode where atoms in- 
side the same layer are "optical" with respect to each 
other. Of course "primed" optical modes do not exist for 
graphene, since there is only one layer (two atoms) per 
unit cell. 

We observe that stacking has a negligible effect on all 
the frequencies above 400 cm~^ , since both rhombohedral 
graphite and hexagonal graphite show nearly the same 
dispersions except for the F-A branch and the in-plane 
dispersions near F. The in-plane part of the dispersions is 
also very similar to that of graphene, except of course for 
the low optical branches (below 400 cm"-'^) that appear 
in graphite and are not present in graphene. 

For graphite as well as diamond GGA tends to make 
the high optical modes weaker while LDA makes them 
stronger than experimental values. The opposite hap- 
pens for the low optical modes, and for the F-A branch 
of graphite; the acoustic modes show marginal differences 
and are in very good agreement with experiments. Over- 
all, the agreement of both LDA and GGA calculations 
with experiments is very good and comparable to that 
between different measurements. 

Some characteristic features of both diamond and 
graphite are well reproduced by our ab-initio results, such 
as the LO branch overbending and the associated shift 
of the highest frequencies away from F. Also, in the 
case of graphite, rhombohedral graphite and graphene, 
the quadratic dispersion of the in-plane ZA branch in 
the vicinity of F is observed; this is a characteristic fea- 
ture of the phonon dispersions of layered crystals^^'Si, 
observed experimentally e.g. with neutron scattering^. 
Nevertheless, some discrepancies are found in graphite. 
The most obvious one is along the F-M TA branch, where 
EELS^^ data show much higher frequencies than calcula- 
tions. Additionally several EELS experiments^SiS report 
a gap between the ZA and ZO branches at K while these 
cross each other in all the calculations. In these cases the 
disagreement could come either from a failure of DFT 
within the approximations used or from imperfections in 
the crystals used in the experiments. 

There are also discrepancies between experimental 
data, in particular in graphite for the LA branch around 
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TABLE II: 


Structural and elastic properties of 


graphite according 


to LDA, GGA, and experiments 






LDA fully 


GGA fully 


GGA using 


GGA with 


Experiment 




theoretical 


tlicorctical 


cxp. Co 


2^^"^ derivatives 


(300 K) 








in Eqs. il6al 


taken at exp. co/ao 


Lattice constant ao(a.u.) 


4.61 


4.65 


4.65 


4.65(fixed) 


4.65±0.003" 


— ratio 


2.74 


3.45 


3.45 


2.725(fixed) 


2.725±0.001'' 


Cii + Ci2 (GPa) 


1283 


976 


1235 


1230 


1240±40'' 


C33 (GPa) 


29 


2.4 


1.9 


45 


36.5±1 ' 


Ci3 (GPa) 


-2.8 


-0.46 


-0.46 


-4.6 


15±5 ' 


Bo (GPa) 


27.8 


2.4 


1.9 


41.2 


35.8 ° 


C' (GPa) 


225 


164 


207 


223 


208.8 " 



"Refs. 51 52 53, as reported by Ref.|22 
''Ref.ia 

'^Ref.lM as reported by R.ef.l22 
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FIG. 3; GGA ab-initio phonon dispersions (solid lines) and vibrational density of states (VDOS) for diamond. Experimental 
neutron scattering data from Ref. (3 are shown for comparison (circles) . 



TABLE III: Phonon frequencies of diamond at the high- 
symmetry points r, X and L, in cm~^. 





To 


Xta 


Xto 


Xlo 


Lta 


Lla 


Lto 


Llo 


LDA" 


1324 


800 


1094 


1228 


561 


1080 


1231 


1275 


GGA' 


1289 


783 


1057 


1192 


548 


1040 


1193 


1246 


Exp." 


1332 


807 


1072 


1184 


550 


1029 


1206 


1234 



°R.ef..lL 

''Present calculation 
=Ref.|g 



K: EELS data from Ref. ^ agree with our ab-initio re- 
sults while those from Ref. 57 deviate from them. 

Finally, we should stress again the dependence of the 
graphite phonon frequencies on the in-plane lattice pa- 
rameter and c/a ratio. The results we have analyzed 
so far were obtained using the theoretical in-plane lat- 
tice parameter a and the experimental c/a ratio for both 
GGA and LDA. Since the LDA theoretical c/a is very 
close to the experimental one (2.74 vs. 2.725) and the 
interlayer bonding is very weak, these differences do no 
matter. However this is not the case for GGA, as the the- 
oretical c/a ratio is very different from the experimental 
one (3.45 vs. 2.725). Fig. [7| and the second column of 



Table HVl show results of GGA calculations performed at 
the theoretical c/a. Low frequencies (below 150 cm~^) 
between F and A are strongly underestimated, as are 
the ZO' modes between F and M, while the remaining 
branches are barely affected. 

The high-frequency optical modes are instead strongly 
dependent on the in-plane lattice constant. The differ- 
ence between the values of a in LDA and GGA explains 
much of the discrepancy between the LDA optical modes 
and the GGA ones. Indeed, a LDA calculation performed 
at a = 4.65 a.u. and c/a = 2.725 (not shown here) 
brings the phonon frequencies of these modes very close 
to the GGA ones obtained with the same parameters, 
while lower-energy modes (below 1000 cm~^) are hardly 
affected. 

Our final choice to use the theoretical in-plane lattice 
parameter and the experimental c/a seems to strike a 
balance between the need of theoretical consistency and 
that of accuracy. Therefore, the remaining of this section 
is based on calculations performed using the parameters 
discussed above (a = 4.61 for LDA, a — 4.65 for GGA 
and c/a — 2.725 in each case). 

Elastic constants can be extracted from the data on 
sound velocities. Indeed, the latters are the slopes of 
the dispersion curves in the vicinity of F and can be ex- 
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FIG. 4: GGA (solid lines) and LDA (dashed line) ab-initio phonon dispersions for graphite, together with the GGA vibrational 
density of states (VDOS). The inset shows an enlargement of the low-frequency F-A region. The experimental data are EELS 
(Electron Energy Loss Spectroscopy) from Refs|^5, 56, 57 (respectively squares, diamonds , an d filled circles), neutron scattering 
from Ref. Is^ (open circles), and x-ray scattering from Ref. ^3 (triangles). Data for Refs. Issi andl57l were taken from Ref. 
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VDOS 



FIG. 5; GGA ab-initio phonon dispersions for graphene (solid 
lines). Experimental data for graphite are also shown, as in 
Fig. El 



FIG. 6; GGA ab-initio phonon dispersions for rhombohe- 
dral graphite. The inset shows an enlargement of the low- 
frequency F-A region. 



pressed as the square root of linear combinations of elas- 
tic constants (depending on the branch considered) over 
the density (see Ref. |63 for details). We note in pass- 
ing that we computed the density consistently with the 
geometry used in the calculations (see Table Hvl for de- 
tails, first column for LDA and third one for GGA), and 
not the experimental density. Our results are shown in 



Table |V| 

The overall agreement with experiment is good to very 
good. LDA leads to larger elastic constants, as expected 
from the general tendency to "overbind" , but still agrees 
well with experiment. For diamond, the agreement is 
particularly good. As for C13 in graphite, it is quite 
difficult to obtain it from the dispersion curves since it 
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TABLE IV: Phonon frequencies of graphite and derivatives at the high-symmetry points A, F, M and K, in cm ^. The lattice 
constants used in the calculations are also shown. 







Graphite 




Rhombo. graphite 


Graphene 


Graphite 


Functional 


LDA 


GGA 


GGA 


GGA 


GGA 


Experiment 


In-planc lattice ct. ao 


4.61 a.u. 


4.65 a.u. 


4.65 a.u. 


4.65 a.u. 


4.65 a.u. 


4.65 a.u. 


Intcilaycr distance/ ao 


1.36 


1.725 


1.36 


1.36 


15 


1.36 


At A/TO' 


31 


6 


29 






35" 


Ala/lo' 


80 


20 


96 






89" 


Alo 


897 


880 


878 








Ato 


1598 


1561 


1564 








Tlo' 


44 


8 


41 


35 




49° 


^ zo' 


113 


28 


135 


117 




yo , iz6 


^zo 


899 


881 


879 


879 


881 


861 


^lo/to 


1593 


1561 


1559 


1559 


1554 


1590', 1575^ 


1604 


1561 


1567 








MzA 


478 


471 


477 


479 


471 


471°, 465 , 451 


Mta 


630 


626 


626 


626 


626 


630'' 


Mzo 


637 


634 


634 


635 


635 


670' 


Ml a 


1349 


1331 


1330 


1330 


1328 


1290^ 


Mlo 


1368 


1346 


1342 


1344 


1340 


1321 = 


Mto 


1430 


1397 


1394 


1394 


1390 


1388', 1389' 


Kza 


540 


534 


540 


535 


535 


482°, 517°, 530° 


Kzo 


544 


534 


542 


539 


535 


588'*, 627" 


Kta 


1009 


999 


998 


998 


997 




Kla/lo 


1239 


1218 


1216 


1216 


1213 


1184', 1202° 


Kto 


1359 


1308 


1319" 


1319 


1288' 


1313'', 1291° 



''Ref.^ 
'^Ref. 
-^Ref. 
'^Ref. 
•''Ref.l 

^Note that a direct calculation of this mode with DFPT (instead 
of the Fourier interpolation result given here) leads to a significantly 
lower value in the case of graphite — 1297 cm~^ instead of 1319 
cm~^. This explains much of the discrepancy between the graphite 
and graphene result, since in the latter we used a denser q-points 
mesh. This effect is due to the Kohn anomaly occurring at Kbb. 




F VDOS 



FIG. 7: GGA ab-initio phonon dispersions for graphite at 
the theoretical c/a. The inset shows an enlargement of the 
low- frequency F-A region. 



TABLE V: Elastic constants of diamond and graphite as 
calculated from the phonon dispersions, in GPa. 







Diamond 




Graphite 


Functional 


GGA 


Exp. 


LDA 


GGA Exp. 


Cii 


1060 


1076.4 ± 0.2' 


1118 


1079 1060 ± 20" 


Cl2 


125 


125.2 ± 2.3' 


235 


217 180 ± 20" 


C44 


562 


577.4 ± 1.4' 


4.5 


3.9 4.5 ± 0.5" 


C33 






29.5 


42.2 36.5±1 " 



"Ref . M 
''Ref.Eo 



and thus several thermodynamical quantities. Before ex- 
ploring this in Section Hvl we want to discuss the nature 
and decay of the interatomic force constants in carbon 
based materials. 



Interatomic force constants 



enters the sound velocities only in a linear combination 
involving other elastic constants, for which the error is 
almost comparable to the magnitude of C13 itself. 

An accurate description of the phonon dispersions al- 
low us to predict the low-energy structural excitations 



As explained in Section lll Al the interatomic force con- 
stants Ci, j (R — R') are obtained in our calculations from 
the Fourier transform of the dynamical matrix Di^j{q) 
calculated on a regular mesh inside the Brillouin zone 
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FIG. 8: Decay of the norm of the interatomic force constants 
as a function of distance for diamond (thin solid line) and 
graphene (thick solid line), in a semi-logarithmic scale. The 
dotted and dashed lines show the decay for diamond along 
the (100) and (110) directions. 

(8 X 8 X 8 for diamond, 8 x 8 x 4 for graphite and 16x16x1 
for graphene). This procedure is exactly equivalent (but 
much more efficient) than calculating the interatomic 
force constants with frozen phonons (up to 47 neighbors 
in diamond and 74 in graphene) . At a given R, Ci^ j (R) 
is actually a 2"'^ order tensor, and the decay of its norm 
(defined as the square root of the sum of the squares of 
all the matrix elements) with distance is a good measure 
to know the effect of distant neighbors. In Fig. |S| we 
have plotted the natural logarithm of such a norm with 
respect to the distance from a given atom, for diamond 
and graphene. The norm has been averaged on all the 
neighbors located at the same distance before taking the 
logarithm. 

The force-constants decay in graphene is slower than 
in diamond, and it depends much less on direction. In 
diamond decay along (110) is much slower than in other 
directions due to long-range elastic effects along the cova- 
lent bonds. This long-range decay is also responsible for 
the flattening of the phonon dispersions in zincblcnde and 
diamond semiconductors along the K-X line (see Fig. O 
and Ref . 0, for instance) . 

In Fig. 1^ we show the decay plot for graphite and 
graphene, averaged over all directions. The graphite in- 
teratomic force constants include values corresponding to 
graphene (in-plane nearest neighbors) and smaller values 
corresponding to the weak interlayer interactions. 

It is interesting to assess the effects of the truncation 
of these interatomic force constants on the phonon dis- 
persion curves. This can be done by replacing the force 
constants corresponding to distant neighbors by zero. In 
this way the relevance of short-range and long-range con- 
tributions can be examined. The former are relevant for 
short-range force-constant models such as the VFF (Va- 
lence Force Field)^ or the 4NNFC (4*'* Nearest-Neighbor 
Force Constant used e.g. in graphene. Note however 
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FIG. 9: Decay of the norm of the interatomic force constants 
as a function of distance for graphite (thin solid line) and 
graphene (thick solid line). 



that a simple truncation is not comparable to the VFF 
or 4NNFC models, where effective interatomic force con- 
stants would be renormalized. 

Figs. 1101 and 1111 show the change in frequency for se- 
lected modes in diamond and graphene as a function of 
the truncation range. The modes we chose are those most 
strongly affected by the number of neighbors included. 

For diamond, our whole supercell contains up to 47 
neighbors, and the graph shows only the region up to 
20 neighbors included, since the selected modes do not 
vary by more than lcm~^ after that. With 5 neigh- 
bors, phonon frequencies are already near their converged 
value, being off by at worst 4% off; very good accuracy 
(5cm~^) is obtained with 13 neighbors. 

For graphene, our 16x16x1 supercell contains up to 74 
neighbors, but after the 30*'* no relevant changes occur. 
At least 4 neighbors are needed for the optical modes to 
be converged within 5-8%. Some acoustic modes require 
more neighbors, as also pointed out in Ref.|3 As can be 
seen in Fig. 1111 the frequency of some ZA modes in the 
F-M branch (at about one fourth of the branch) oscillates 
strongly with the number of neighbors included, and can 
even become imaginary when less than 13 are used, re- 
sulting into an instability of the crystal. This behavior 
does not appear in diamond. Also, the Kto mode keeps 
varying in going down from 20 to 30 neighbors, though 
this effect remains small (8 — 9cm^^). This drift could 
signal the presence of a Kohn anomaly^''. Indeed, at the 
K point of the Brillouin zone the electronic band gap 
vanishes in graphene, so that a singularity arises in the 
highest optical phonon mode. Therefore a finer q-point 
mesh is needed around this point, and longer-ranged in- 
teratomic force constants. This effect is discussed in de- 
tail in Ref. ,2^ 
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FIG. 10: Phonon frequencies of diamond as a function of the 
number of neighbors included in the interatomic force con- 
stants: To (solid line), Xto (dotted line), and Lta (dashed 
line). 
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FIG. 11: Phonon frequencies of graphene as a function of 
the number of neighbors included in the interatomic force 
constants: Tj^o/to (solid line), Kto (dot-dashed), Mzo 
(dashed), and for the dotted line a phonon mode in the ZA 
branch one-fourth along the F to M line. 



IV. THERMODYNAMICAL PROPERTIES 

We present in this final section our results on the 
thermodynamical properties of diamond, graphite and 
graphene using the quasi-harmonic approximation and 
phonon dispersions at the GGA level. As outlined in Sec- 
tion we first perform a direct minimization over the 
lattice parameter(s) {a.i} of the vibrational free energy 
F{{ai}, T) (Eq.lHI). This gives us, for any temperature T, 
the equilibrium lattice parameter (s), shown in Figs. IT^ 
1131 and [T^ For diamond and graphene, we used in Eq. [S] 
the equations of state obtained from the ground state 
calculations presented in Section IlII Al For graphite this 
choice would not be useful or accurate, since the theoret- 
ical c/a is much larger than the experimental one. So we 
forced the equation of state to be a minimum for a=4.65 
a.u. and -=2.725 (fixing only c/a and relaxing a would 
give a=4.66 a.u., with negligible effects on the thermal 




Temperature (K) 

FIG. 12: Lattice parameter of diamond as a function of tem- 
perature 



expansion). In particular, our "corrected" equation of 
state is obtained by fitting with a fourth order polyno- 
mial the true equation of state around the experimental 
a and c/a, and then dropping from this polynomial the 
linear order terms. Since the second derivatives of the 
polynomial are unchanged, this is to say we keep the 
elastic constants unchanged. The only input from exper- 
iments remains the c/a ratio. We have also checked the 
effect of imposing to C13 its experimental value (C13 is 
the elastic constant that is less accurately predicted), but 
the changes were small. 

The dependence of the phonon frequencies on the lat- 
tice parameters was determined by calculating the whole 
phonon dispersions at several values and interpolating 
these in between. For diamond and graphene we used 
four different values of a (from 6.76 to 6.85 a.u. for di- 
amond, and from 4.654 to 4.668 a.u. for graphene) and 
interpolated them with a cubic polynomial. For graphite, 
since the minimization space is two dimensional, we re- 
stricted ourselves to a linear interpolation and calculated 
the phonon dispersions at three different combinations of 
the lattice constants : (a, c/a) = (4.659,2.725), (4.659,2.9) 
and (4.667,2.725). 

Before considering thermal expansion, we examine the 
zero-point motion. Indeed, lattice parameters at K 
are different from their ground state values. The effects 
of the thermal expansion (or contraction) up to about 
1000 K are small compared to the zero-point expansion 
of the lattice parameters. In diamond, a expands from 
6.743 a.u. (ground state value) to 6.768 a.u., a difference 
of 0.4%. For graphene, a is 4.654 a.u. at the ground 
state and 4.668 a.u. with zero-point motion corrections 
(-1-0.3%); for graphite a increases from 4.65 to 4.664 a.u. 
(-hO.3%) and c from 12.671 to 12.711 (-hO.3%). The in- 
crease is similar in each case, and comparable to the dis- 
crepancy between experiments and GGA or LDA ground 
states. 

The coefficients of linear thermal expansion at any T 
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FIG. 13; In-plane lattice parameter of graphite (solid line) 
and graphene (dashed line) as a function of temperature 
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FIG. 14: Out-of-plane lattice parameter of graphite as a func- 
tion of temperature 



are obtained by numerical differentiation of the previous 
data. Results are shown in Figs . IT^ IT^ and [T7I 

For the case of diamond, we have also plotted the lin- 
ear thermal expansion coefficient calculated using the 
Griineisen formalism (Eq. I12|l instead of directly mini- 
mizing the free energy. While at low temperature the 
two curves agree, a discrepancy becomes notable above 
1000 K, and direct minimization should be performed. 
This difference between Griineisen theory and direct min- 
imization seems to explain much of the discrepancy be- 
tween the calculations of Ref. JT, and our results. Fi- 
nally a Monte-Carlo path integral study by Herrero and 
RamireaS^, which does not use the QHA, gives very sim- 
ilar results. 

For graphite, the in-plane coefficient of linear thermal 
expansion slightly overestimates the experimental values, 
but overall the agreement remains excellent, even at high 
temperatures. Out-of-plane, the agreement holds well up 
to 150 K, after which the coefficient of linear thermal 



FIG. 15: Coefficient of linear thermal expansion for diamond 
as a function of temperature. We compare our QHA-GGA ab- 
initio calculations (solid line) to experiments (Ref. 10, filled 
circles), a path integral Monte-Carlo study using a Tersoff 
empirical potential (Ref. 1651 open squares) and the QHA- 
LDA study by Pavone et aii^ (dashed line). The QHA-GGA 
thermal expansion calculated using the Griineisen equation 
fEg. I12II is also shown (dotted line). 
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FIG. 16: In-plane coefficient of linear thermal expansion as a 
function of temperature for graphite (solid line) and graphene 
(dashed line) from our QHA-GGA ab-initio study. The ex- 
perimental results for graphite are from Ref. El (filled circles) 
and Ref. (open diamonds) . 



expansion is underestimated by about 30% at 1000 K. 

In-plane, the coefffcient of linear thermal expansion is 
confirmed to be negative from to about 600 K. This 
feature, absent in diamond, is much more apparent in 
graphene, where the coefficient of linear thermal expan- 
sion keeps being negative up to 2300 K. This thermal con- 
traction will likely appear also in single- walled nanotubes 
(one graphene sheet rolled on itself)^. Some molecular 
dynamics calculations^itSi have already pointed out this 
characteristic of SWNT. 

To further analyze thermal contraction, we plotted 
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FIG. 17; Out-of-plane coefficient of linear tfiermal expansion 
as a function of temperature for graphite from our QHA-GGA 
ab-initio study (solid line) . The experimental results are from 
Ref. [3 (filled circles) , and Ref. Q (open diamonds) . 
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FIG. 18: Ab-initio mode Griineisen parameters for diamond. 

in Figs. UHI HOI and 1^ the mode Griineisen pa- 
rameters (see Section III Bp of diamond, graphene and 
graphite. These have been obtained from an interpo- 
lation of the phonon frequencies by a quadratic (or hn- 
ear, for graphite) polynomial of the lattice constants, and 
computed at the ground state lattice parameter. 

The diamond Griineisen parameters have been already 
calculated with LDA (see Refs. HMO); our GGA re- 
sults agree very well with these. In particular, all 
the Griineisen parameters are shown to be positive (at 
odds with other group IV semiconductors such as Si 
or Ge). The situation is very different in graphite and 
graphene, where some bands display large and nega- 
tive Griineisen parameters (we have used the definition 

-v. (a) - ° du>jiq) ^ 

lA'i) — 2wj(q) da >■ 

While not visible in the figure, the Griineisen parame- 
ter for the lowest acoustic branch of graphite becomes as 
low as -40, and as low as -80 for graphene. Therefore, at 
low temperatures (where most optical modes with posi- 




FIG. 19: Ab-initio in-plane mode Griineisen parameters for 
graphite. 
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FIG. 20: Ab-initio mode Griineisen parameters for graphene. 

five Griineisen parameters are still not excited) the con- 
tribution from the negative Griineisen parameters will be 
predominant and thermal expansion (from Eq. I12|l nega- 
tive. 

The negative Griineisen parameters correspond to the 
lowest transversal acoustic (ZA) modes, and in the case 
of graphite to the (ZO') modes, which can be described 
as "acoustic" inside the layer and optical out-of-plane 
(see Section IIII Bll . Indeed, the phonon frequencies for 
such modes increase when the in-plane lattice parameter 
is increased, contrary to the usual behavior, because the 
layer is more "stretched" when a is increased, and atoms 
in that layer will be less free to move in the z direction 
(just like a rope that is stretched will have vibrations of 
smaller amplitude, and higher frequency). In graphite 
these parameters are less negative because of the inter- 
action between layers: atoms are less free to move in the 
z-direction than in the case of graphene. 

This effect, known as the "membrane effect" , was pre- 
dicted by Lifshitz^^ in 1952, when he pointed out the 
role of these ZA modes (also called "bending modes") 
in layered materials. In particular, several recent stud- 
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FIG. 21: Ab-initio out-of-plane mode Griineisen parameters 
for graphite. 



ies have highlighted the relevance of these modes to the 
thermal properties of layered crystals such as graphite, 
boron nitride and gallium sulfideS2iS242£. 

The knowledge of the equilibrium lattice constant (s) 
at any temperature allows us also to calculate the de- 
pendence of elastic constants on temperature. To do so 
we calculated the second derivatives of the free energy 
(Eq. |HJ| vs. lattice constant(s) at the finite-temperature 
equilibrium lattice parameter (s). We checked that this 
was equivalent to a best fit of the free energy at T around 
the equilibrium lattice parameter(s). 

Results are shown in Figs. |221 and [23 (diamond and 
graphite respectively). Again, the zero-point motion has 
a significant impact on the elastic constants; the agree- 
ment with experimental data for the temperature depen- 
dence of the ratio of the bulk modulus of diamond to its 
298 K value is excellent (upper panel of Fig. I22|) . 

We note that the temperature dependence of the bulk 
modulus of diamond has already been obtained by Karch 
et af^ using LDA calculations. 

As final thermodynamic quantities, we present results 
on the heat capacities for all the systems considered, at 
constant volume (C^) and constant pressure (Cp). Cy 
has been computed using Eq. ^1 in which we used at 
each temperature T the interpolated phonon frequencies 
calculated at the lattice constant (s) that minimize the 
respective free energy. To calculate Cp, we added to C„ 
the additional term Cp — = TVoBoay where Vq is the 
unit cell volume, ay the volumetric thermal expansion 
and Bq the bulk modulus. All these quantities were taken 
from our ab-initio results and evaluated at each of the 
temperatures considered. The difference between Cp and 
Cv is very small, at most about 2% of the value of C^ 
for graphite and 5% for diamond. Note that Cp and Cy 
shown on the figures are normalized by dividing by the 
unit cell mass. 

The heat capacity of diamond, graphite and graphene 
are almost identical except at very low temperatures. 
Agreement with experimental data is very good. 



FIG. 22; Lower panel: Bulk modulus Bo{T) of diamond as a 
function of temperature. The filled circle indicates the value 
of the bulk modulus (as in TablePi before accounting for zero- 
point motion. Upper panel: theoretical (solid line) and ex- 
perimental values (Ref.^2j open circles) for the ratio between 
Bo{T) and _Bo(298-ft') in the low temperature regime. 
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FIG. 23: Elastic constants of graphite (Cn -)- C12, C13, C33) 
and bulk modulus (-Bo) as a function of temperature. The 
filled circles (at K) indicate their ground state values (as in 
Table ^ before accounting for zero-point motion. 



V. CONCLUSION 

We have presented a full ab-initio study of the struc- 
tural, vibrational and thermodynamical properties of di- 
amond, graphite and graphene, at the DFT-GGA level 
and using the quasi-harmonic approximation to derive 
thermodynamic quantities. All our results are in very 
good agreement with experimental data: the phonon dis- 
persions are well-reproduced, as well as most of the elastic 
constants. In graphite, the C33 elastic constant and the 
F to A phonon dispersions (calculated here with GGA for 
the first time) were found to be in good agreement with 
experimental results provided the calculations were per- 
formed at the experimental c/a. Only the C13 constant 
remains in poor agreement with experimental data. 

The decay of the long-ranged interatomic force con- 
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FIG. 24; Constant pressure heat capacity for diamond (solid 
line). Experimental results are from Refs . l49l and (circles), 
as reported by Ref. IbrL 
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FIG. 25: Constant pressure heat capacity for graphite (solid 
line). Experimental results are from Ref. (squares), as 
reported by Ref. 75. 
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stants was analyzed in detail. It was shown that interac- 
tions in the (110) direction in diamond are longer-ranged 
than these in other directions, as is characteristic of the 
zincblende and diamond structures. For graphene and 
graphite, in-plane interactions are even longer-ranged 
and phonon frequencies sensitive to the truncation of the 
interatomic force constants. 



Thermodynamical properties such as the thermal ex- 
pansion, temperature dependence of elastic moduli and 
specific heat were calculated in the quasi-harmonic ap- 
proximation. These quantities were all found to be in 
close agreement with experiments, except for the out- 
of-plane thermal expansion of graphite at temperatures 
higher than 150 K. Graphite shows a distinctive in- 
plane negative thermal-expansion coefficient that reaches 
the minimum around room temperature, in very good 
agreement with experiments. This effect is found to be 
three times as large in graphene. In both cases, the 
mode Griineisen parameters show that the ZA "bending" 
acoustic modes are responsible for the contraction, in a 
direct manifestation of the membrane effect predicted by 
LifshitaSi in 1952. These distinctive features will likely 
affect the thermodynamical properties of single-walled 
and multiwall carbon nanotubes^i*2SiSi. 
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FIG. 26: Constant volume heat capacity for graphite (solid 
line) , graphene (dashed line) and diamond (dotted line) . The 
inset shows an enlargement of the low temperature region. 
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